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We have computed the complete next-to-next-to-leading order (NNLO) contributions to the splitting functions 
governing the evolution of unpolarized parton densities in perturbative QCD. Our results agree with all partial 
results available in the literature. We illustrate the methods used for this calculation with some examples and 
display selected results to show the size of the NNLO corrections and their effect on the evolution. 



1. INTRODUCTION 

Parton distributions form indispensable ingre- 
dients for the analysis of all hard-scattering pro- 
cesses involving initial-state hadrons. The scale- 
dcpcndcncc of these distributions can be de- 
rived from first principles in terms of an ex- 
pansion in powers of the strong coupling con- 
stant a s . The corresponding nth-order co- 
efficients governing the evolution are referred 
to as the n-loop splitting functions. Includ- 
ing the terms up to order in the evo- 
lution of parton distributions, together with 
the corresponding results for the hard partonic 
cross sections of a given observable, one obtains 
the N™LO (leading-order, next-to-lcading-order, 
ncxt-to-next-to-leading-order, etc.) approxima- 
tion of perturbative QCD. 

The standard approximation for most impor- 
tant processes is presently the next-to-leading or- 
der, the corresponding one- and two-loop split- 
ting functions being known for a long time. How- 
ever, the NNLO corrections need to be included 
in order to arrive at quantitatively reliable pre- 
dictions for hard processes at present and future 
high-energy colliders. 

In a series of recent papers |1I2I3| , we have pub- 
lished the complete unpolarized three-loop split- 
ting functions. Here, we present some aspects of 
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the methods and selected results. 

2. METHODS 

The method of calculation employs the opti- 
cal theorem, which relates the total cross section 
for a given process to the imaginary part of the 
forward Compton amplitude. We study deep- 
inelastic scattering of a boson with Euclidean (off- 
shell) momentum Q, off a parton with (on-shell) 
momentum P and apply the operator product 
expansion. In the Bjorkcn limit, Q 2 — > oo and 
x = Q 2 /(2P ■ Q) fixed, this yields a relation 
between matrix elements of parton operators of 
leading twist and Mellin moments of Feynman 
diagrams contributing to the forward Compton 
amplitude. Specifically, using dimensional regu- 
larization in D = 4 — 2e, we are able to com- 
pute the anomalous dimensions 7 of the parton 
operators, i.e. the integer-A Mellin moments of 
splitting functions P, 



7 (A) 



dx x N ~ 1 P{x) 



(1) 



from the divergence in e of the A-th Mellin mo- 
ment of the corresponding Feynman diagrams. 

In general, the (anti-)quark (anti-)quark split- 
ting functions, constrained by charge conjugation 
invariance and flavour symmetry, are given for 
flavours i, k by 
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They can be composed into three independent 
types of non-singlet splitting functions, 



P ± 

ns 

and 



pv ± pv 

qq qq 



(3) 
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which govern the evolution of the quark flavour 
asymmetries q^ s and the valence distribution q^ s , 
respectively, 

d 
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±,V. 
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The quark-quark singlet splitting function P qq 
is expressed as 



P 



qq 



^nt+^/OPqq + ^qq) = ^ 
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and gluon-quark and quark-gluon splitting func- 
tions arc given by 



P qg — 71/ P qig 



(7) 



in terms of the flavour-independent splitting func- 



tions P n 



p 

singlet quark distribution q s and the gluon distri- 
bution g evolve according to 



and P„, 
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Eqs. (JSJ and (jSJ form the well-known 2n/ — 1 
scalar non-singlet and 2x2 singlet evolution equa- 
tions, where (g> stands for the Mcllin convolution. 

In the present calculation, for the complete set 
of singlet and non-singlet splitting functions, we 
had to compute 18 one-loop Fcynman diagrams at 
LO, 350 two-loop diagrams at NLO, and a total of 
9607 three- loop diagrams at NNLO. All diagrams 
were generated automatically with the diagram 
generator Qgraf and for all symbolic manip- 
ulations we used the latest version of FORM 0. 
The NNLO calculation required significant en- 
hancements of the capabilities of Form [Hj . 

In order to illustrate the details of the calcu- 
lation, let us pick a particular Fcynman diagram 
at three loops, which contributes to the forward 
Compton amplitude for the scattering of a photon 
off a quark. 




Here and below, the fat line indicates the flow of 
the quark momentum P through the diagram. 

In a first step, we need a classification of the 
topology of the loop integral. This is done in two 
steps. Any loop integral with external momenta 
P and Q, we classify according to 

1. the topology of the underlying "two-point 
function" if we nullify P. 

2. the flow of the parton momentum P. 

The first step yields at the top-level the topology 
types ladder, benz and non-planar, where we use 
the notation of Refs. |7I8| . The second step distin- 
guishes, within each topology, between so-called 
basic building blocks (BBB) which have a simple 
P-flow with one P-dependent propagator only, 
and so-called composite building blocks (CBB), 
which have a more complicated P-flow. At the 
top-level, there are 10 BBB (3 ladder, 5 benz and 
2 non-planar) and 32 CBB (10 ladder, 16 benz 
and 6 non-planar). The smaller number of non- 
planar topologies is due to symmetries. For illus- 
tration purposes, we represent the complete set 
of top-level BBB by pictograms as follows. 




A quick inspection shows, that our example is 
a particular ladder type BBB, 



3 



(2P-Q) 



N 



The equation indicates that we calculate the 
7V-th Mellin moment of this diagram. This is 
precisely the dimensionless coefficient Cat given 
on the right hand side. On the left hand side, the 
fat line in the pictogram represents the flow of the 
momentum P, in form of only one internal prop- 
agator containing both, a loop momentum, say 
li, and the external quark momentum P. Thus, 
on the left hand side the N-th term in a Taylor 
expansion generates the contribution to the N-th 
Mellin moment, 



(2P • hf 



(2P-h) 



N 



2\N 



(9) 



Thus far, the set-up is completely analogous to 
the calculation of the lowest six/seven (even or 
odd) integer- N Mellin moments of the three-loop 
splitting functions [911011 lj . where the Mincer 
program |7I8| was used as a tool to solve the inte- 
grals. As a new feature here we are dealing with 
symbolic N. This was tested before up to two 
loops 22- At three loops, it leads, for instance, 
to integrals of the type 

(2P-h) k 1 



Jd D l, 



i 2 

2 



where the powers n, k are symbolic. 

However, we can switch to non-symbolic (fixed) 
positive powers n, k and values of N at any point 
of the derivations and calculations, after which 
the Mincer program can be invoked to verify 
that the results are correct. From a practical 
point of view this is the most powerful feature 
of the Mellin-space approach, as it allows for ex- 
tremely efficient checks. 

Having classified the integrals, we actually need 
to solve them, too. To this end, let us start the 
discussion with the most general form of a three- 
loop integral in one of the top-level topologies lad- 
der, benz or non-planar, 



I(N, ft, V, k) = 
(2P-h) Kl 



n 



(10) 



(2P-l a ) K *(2l g -l 10 ) 



(l 2 )^((P-h) 2 )^ . . . (Z|)M8 ((P - Z 8 )2)*s 



Here, fx = Hi, . . . , fj,8, v = v\,...,v% and k = 
K\, . . . , Kg are symbolic parameters. Further- 
more, the li are for i = 1,2,3 independent loop 
momenta. For i = 4, . . . , 8 they can be ex- 
pressed by the former and the external momen- 
tum Q, the precise relations depending on the 
topology. The momenta lg,ho denote an irre- 
ducible scalar product between the momenta Q 
and/or ij ,i = 1, . . . , 8, again the precise relation 
being topology depended. 

Applying relations based on integration by 
parts [13114115116] . scaling equations, form- factor 
analysis ^Zj an d some equations that fall in a 
special category because they involve higher twist 
and a careful study of the parton-momcntum 
limit P ■ P — > 0, we arrive at a system of linear 
equations for a given integral I(N, ft, V, k) under 
consideration. 

Solving this linear system amounts to finding 
a scheme in which I(N, ft, v, K) is mapped to a 
set of master integrals and integrals of simpler 
topologies. The general strategy applies two basic 
rules of mapping, 



1. non-planar — > benz 



ladder. 



2. CBB 



BBB. 



Here, the first rule is understood to hold for com- 
mon sub-topologies. For instance, the non-planar 
topology contains common sub-topologies with 
the benz and the ladder topology [718] , 

In an operator approach, this can be realized 
by diagonalizing the linear system for symbolic 
indices. One thus obtains lowering operators for 
individual v i or /Xj , as well as recursion relations in 
the Mellin moment N. The latter constitute dif- 
ference equations, which may generally be written 
as 



a (N) I(N) + ai(N) I(N - 1) 
+ ... + a m (N)I{N-m) 



(11) 



G(N) 



To illustrate the latter, let us give an extremely 
simple example for a single-step difference equa- 
tion in TV, which occurs in the reduction of a par- 
ticular type of CBB, 



<n>- 



N+2 



iV+3+3e 2P-Q 
N+2 Q 2 



In the pictogram, all lines denote propagators of 
unit power, except the one with a blob, which 
has to be taken to second power. Again, it is 
understood that the equation holds for the N-th 
Mcllin moment of the diagrams and the fat lines 
indicate propagators with the momentum P. 

Employing the notation of Eq. (|11J) . we can 
write the single-step difference equations as 

G(«).(12) 



N+2 



N+2 



As a remark on the side, imagine for a moment, 
the function G(N) on right hand side would be 
multiplied by an additional factor e _1 . If present, 
such a so-called spurious pole in e would make 
Eq. 1|12J) useless. It would ruin the accuracy of 
the expansion when working only a to given cut- 
off in powers of e. (We do so both for reasons of 
economy and because we cannot evaluate some in- 
tegrals easily beyond certain powers in e.) Thus, 
spurious poles have to be avoided and one of the 
greatest difficulties in deriving reduction schemes 
is to indeed avoid them. 

Eq. 1|12|) can be solved in closed form, 



z W = (-i)* n ^ +3+3£) 



7(0) 



(13) 



N 



t (-D"E (-d j n5 b b ' +3 : 3c w 



i=l 



nf=i(i+2) 



Eq. (|13fl is an example for the occurrence of 
nested sums in the calculation. Its solution re- 
quires as input the boundary value 1(0), which 
can be obtained with the Mincer program. The 
inhomogeneous term G(N) is assumed to be al- 
ready known. It has to be calculated by similar 
means, i.e. reductions and recursions. The suc- 
cessive way of solving difference equations induces 
a strict hierarchy for all topology classes in the re- 
duction scheme. 

As a matter of fact, Eq. (|13|l is a special 
case of a general recursion which was derived 



for I(N, p,, i7, it) with symbolic indices, like most 
other recursions as well. As such, they allow for 
an efficient implementation in Form resulting in 
a largely automatic build-up of nested sums. 

The solution of nested sums as in Eq. I|13|l re- 
sults in harmonic sums |18I19I20I21I22| which are 
recursively defined as 



S± mu ...,m k (N) =i>2^-S' 



m 2 ,...,m k (i) ■ (14) 



In particular, we use four main algorithms for 
harmonic sums, which rely on the underlying al- 
gebra |23I24| . The algorithms express products 
or sums of nested sums again in the basis Eq. l(Tf|) 
of harmonic sums. Specifically, they act on prod- 
ucts, 



(N) S ni ,..., ni (N) , 
sums involving j and N — j 

JV-l 



C ' Ami 



i=i 



3 

x S n 



1 



(N - j)™* 



(15) 



(16) 



conjugations 

N 



N / \ 

E 0) -^s m2 ,..., mk u) 



(17) 



and on sums involving binomials, j and N — j, 

N-l , * 

-E (j) (-iy -^s m2 _ mk (j) (is) 



i 



(N - j)" 1 



S n2 ,..., ni (N - j). 



The solution to Eq. Ijl3(l up to order e _1 as 
required for the calculation of the splitting func- 
tions can be obtained, 

I(N)-( lV vlf4 Si(iV+l) 8 S 1 (N+1) 
I(N)-(-l) N + 1 +3 (N+1)2 ( 19 ) 

4S*i(A + 2) 4 51(^ + 2) 4 



3 N+2 



3 (N+2) 



• gSl(A0 



(m i 2 S 2 (N + 1) , 2 S 2 (N + 2) 
+^S h2 (N) + 3 N+1 + 3 N+2 



5 



-2S 2 (iV) - -NS 2 (N) + 4S 2>1 (N) 
+-NS 2 ,i(N) ~ QS 3 {N) - 2NS 3 (N) 



1 4- 1 



1 



3(iV+l) 2 (N+l) 3 3 (N + 2) 2 
1 \ , 

c 



(7V + 2)= 



N+l 



5i(JV + l) 20 5i(JV + l) Si(iV+2) 

— lo- 



3 (N+l) 2 3 (iV + 1) 3 " iV+2 
44 Si(iV + 2) 10 ^ + 2) 

T (tv+2)2 - T (iv+2)3 - 20Sl{N > 

8S u (iV+l) 8 5i,i(AT+l) 8 5i,i(AT+2) 



3 N+l 



3 (N + iy 



3 N+2 



10 5i, 2 (JV+2) 



- 165i, 2 (JV) - 4JVSi, 2 (JV) 



3 7V+2 
+145'i,2,i(iV) + 4iVSi, 2 ,i(A0 - 24S 1 , 3 (7V) 

rat* (An 58 S 2 (iV + l) 40^V + 1) 
-6^,3^) - y^y- - 

46S 2 (iV + 2) S , 2 (A^+2) 56 



3 7V+2 



-6- 



(iV+2) 



S 2 (N) 



-20NS 2 (N) + 10 S ^{ N ^+6 S ^ N+2) 



N+l 



N+2 



-yS2,lW - yiV5 2 ,i(iV) + y5 2jlil (iV) 

8 62 22 

+ gMS 2 ,M(#) - y^2, 2 (iV) - — iVS 2)2 (iV) 

-12^±g) + 765 3 (iV) 
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+ ™NS 3 (N) - 10S 3 ,!(N) - yJVSyfJV) 
+36S 4 (A0 + 12NS 4 (N) + 32 ——. 

[N+iy 



+- 



164 1 



3 (N + l) 3 (N + l) 



2A ,, 1 — + 16 



(N+2f 



+ ■ 



82 1 
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3 (7V + 2) 3 (N+2) 



The result in Eq. I|19(l is neither short nor in- 
expensive to calculate. Moreover, each integral 
is typically used many times. Therefore, to save 
computer time and disk space, we have tabulated 
large numbers of CBB and BBB integrals. Only 



this tabulation rendered the calculation feasible 
with current computing resources. For the com- 
plete project we have collected table-bases with 
more than 100.000 integrals and a total size of 
tables of more than 3 GBytes. 

Subsequently, the database of integrals was 
used for the calculation of all Fcynman diagrams 
yielding the unrenormalizcd results in Mcllin- 
space in terms of the invariants determined by 
the colour group |25| . harmonic sums and the 
values £3, £4, (5 of the Riemann £- function. 
In physics results the terms with £4 cancel in 
iV-space. The renormalization was carried out 
in the MS-scheme |26l27j and the procedure is 
again the same as for the fixed-moment calcula- 

tions mam. 

3. RESULTS 

Now we present the anomalous dimensions 
7(a s , N) up to the third order in the running cou- 
pling constant a s , the expansion coefficients being 
normalized as 



n=0 



4-7T/ 



(20) 



Our analytical results have been presented in 
Refs. 1 11213) and are too long to be reproduced 
here. We agree with all partial results available 
in the literature, in particular we reproduce the 
lowest six even integer singlet and seven even/odd 
integer non-singlet moments 2 computed before 

Ennnu. 

The numerical results for the singlet anoma- 
lous dimension, i.e. the Mellin transforms of the 
matrix entries in Eq. (JHJ are illustrated in Fig. ^ 
In the top row of Fig. ^ wc show the pertur- 
bative expansion of the diagonal anomalous di- 
mensions 7 qq (AT) and 7 gg (iV) for four flavours at 
a s = 0.2. The pure-singlet (ps) contribution to 
7 qq as defined in Eq. © is displayed separately. 
The bottom row of Fig. ^ shows the off-diagonal 
anomalous dimensions 7 qg (7V) and j gq (N). For 
all cases, the NNLO corrections are significantly 
smaller than the NLO contributions and amount 



(2)H 



2 The moment N = 16 of 7 
puted as an additional check 1281 . 



has recently been com- 



6 





Figure 1. The perturbative expansion of the sin- 
glet anomalous dimensions j(a a ,N). 



Figure 2. The three- loop singlet splitting func- 

(2) 

tions P ab with the leading small- a; terms (dot- 
ted) and the fixed- moment estimates (dashed). 



to less than 2% and 1% for the large diagonal 
quantities 7 qq and 7 gg , respectively, for N > 2. 

In Bjorkcn x-spacc, the N ra LO splitting func- 
tions pW(i) in 

P(a s ,x) = J2 (^y +1 P (n) (x), (21) 

can be obtained from Eq. JQ| and expressed in 
terms of harmonic polylogarithms [29130131 ] by 
means of an inverse Mellin transformation. This 
is a completely algebraic procedure based on the 
fact that harmonic sums occur as coefficients of 
the Taylor expansion of harmonic polylogarithms. 

Again, the analytical results have been given 
in Refs. |1I2I3| . They agree with the available re- 
summation predictions pJS 33 34 35] for the lead- 
ing small- a; logarithms, and those for large-n/ re- 
sults |36I37| . In addition, we have also provided 
easy-to-use accurate parameterizations. 

Our results respect the supersymmetric rela- 
tion between all four singlet splitting functions for 



Ca = G'f =71/ to the extend expected for the 
MS scheme. At large x we agree with Refs. [38 39 
and determine for P^' v , P m and P gg the coeffi- 
cients of the leading 1/(1 — x)+ terms. We find 
that the coefficients of the leading integrable term 
ln(l — x) at order n = 2, 3 are proportional to 
the coefficient of the +- distribution 1/(1 — x) + 
at order n — 1, a result that seems to point to a 
yet unexplored structure. Furthermore, we verify 
the expected simple relation between the leading 
1/(1 — x)+ terms of P qq and P gg . 

In Fig. we plot the three-loop singlet split- 
ting functions of Eq. (jHJ in the small- a; region, 
where the leading small- a; behaviour is of the type 
a; -1 In a;. This agrees with the prediction of the 
leading logarithmic BFKL equation [40141142) . In 
the top row of Fig. [21 we give the three- loop split- 

(2) 

ting functions P ps (pure-singlet quark-quark) 
and P qg (gluon-quark) for four flavours, multi- 
plied by x for display purposes. Also shown are 
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Figure 3. The non-singlet three-loop splitting 
functions P n ~ (^-independent part) and P n |/n./. 



the uncertainty band derived in Ref. 03] using 
the lowest six even-integer moments (lUlllj and 
the leading small- a; terms [22 ■ In the bottom row 
of Fig. |21 we show the three-loop splitting func- 



tions Pgq' (quark-gluon) and P gg ^ (gluon-gluon). 

(2) 

For Pgq , the leading small- a; contribution was 

(2) 

unknown before the present calculation, for P gg 
the leading small- x term has been first obtained 
in Ref. [33]. P gg as a diagonal quantity has been 
additionally multiplied by (1 — x). As illustrated 
in Fig. [3 the leading small- a; terms alone do not 
provide good approximations of the full results 
at experimentally relevant small values of x. At 
x = 10 -4 , for example, they exceed the exact val- 

(2) 

ues of P^ b (x) by factors between 1.6 and 2.0 for 
n t = 4. 

In Fig.|3Jwe display the three- loop non-singlet 



(2) 



splitting functions P^\x) and Pnt >s - On the left 
hand side of Fig. |3J the n/-independent three- 

(2) 

loop contribution P_ q (x) to the splitting func- 
tion P n 7(x), multiplied by (1 — x) for display pur- 
poses is shown. Also shown is a comparison of our 
exact result to the small- X expansion in powers of 
In a;. Here, the leading small- a; terms of the type 
In 4 a; were known before for P^ >± {x) [34l35j . As 
can be seen from Fig.|3]on the left, the coefficients 
of ln fe x increase sharply with decreasing power k. 
Including all logarithmically enhanced terms, one 
still underestimates the complete result by a fac- 
tor as large as 2.0 for Pns at x = 10~ 4 . 



(2)3 



(2)s 

The three- loop contribution P,; s exhibits a 
new colour structure d abc d a b c /n c which appears 
for the first time at three loops. It is due to Feyn- 
man graphs of the following type, involving axial 
currents (the quarks couple to VF-bosons), 
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Recall that at one and two loops P^ B and Pnl^ B 
both vanish. PnV s is displayed in Fig. |3J on the 

(2)s 

right. Quiet unexpectedly, Pns also behaves like 
In x for x — > 0, and here the leading small- a; 
terms do indeed provide a reasonable approxima- 
tion. In fact, this function dominates the small- a; 
behaviour of the non-singlet splitting functions, 
for Tij = 4 being, for example, about 7 times 

larger than Pi 2 ' ± (a;) at x = 10~ 4 . 

Let us next illustrate the numerical effect of the 
three-loop splitting functions P^ b ' (x) on the evo- 
lution of the singlet-quark and gluon distributions 
q s {x, and g(x, M/)- F° r ai l figures we choose a 
reference scale [i ? = fx§ ~ 30 GeV 2 - a scale rel- 
evant, for example, for deep-inelastic scattering 
both at fixed-target experiments and the ep col- 
lider HERA - and employ the sufficiently realistic 
model distributions 

xq s (x,^)= (22) 
0.6.T- - 3 (l-a;) 3 - 5 (l + 5.0a; - 8 ), 

x9(x,Vo)= (23) 
1.6.T-°- 3 (l-.T) 4 - 5 (l-0.6a; - 3 ), 

irrespective of the order of the expansion to facil- 
itate the comparison of the LO, NLO and NNLO 
contributions to the splitting functions. 

In Fig. 01 we show in the top row the pertur- 
bativc expansion of the scale derivative q s = 
(ilng s /dln/i| of the singlet quark distribution, 
i.e. the top row of Eq. JHJl, at /i? = fi§ for 
four flavours, a s = 0.2, and the initial conditions 
specified in Eqs. l|25fl. The bottom row of 
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Figure 4. The perturbativc expansion of the 
scale derivatives |jSJ of the singlet distribu- 
tions (221, (mi- 



Fig. 0] shows the same for the gluon distribution 
g, i.e. for the bottom row of Eq. ||SJ). The spikes 
close to x = 0.1 in the right parts of both figures 
are due to zeros of the LO and NLO predictions 
and do not represent large corrections. 

The NNLO corrections are small at large x with 
respect to both the total derivative and the NLO 
contributions. At small- a; the NLO contributions 
are very large for the quark evolution. Conse- 
quently the total NNLO corrections, while reach- 
ing 10% at x = 10 ~ 4 , remain smaller than the 
NLO results by a factor of eight or more over the 
full arrange. 

In Fig.JSlwe show the perturbative expansion of 
the logarithmic scale derivative dlng~ s /<iln/ij 
for a characteristic non-singlet quark distribution 



,0.5 



(1-xf 



±,v, 



(24) 



at the standard scale fx r = fj,f. In addition, on the 
right hand side of Fig. [3] the non-singlet quark 



d In q NS / d In u| 
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Figure 5. The perturbative expansion of the scale 
derivative (jSJ) of the non-singlet distributions g_ 
and q v for the input l|24[) . 



distribution q^ s is also displayed. 

Finally, we turn to the stability of the pertur- 
bative expansions under variations of the renor- 
malization scale fj, r - For fi r ^ [if the expansion 
of the splitting functions in Eq. (|21|l is, using the 
abbreviation a s = a s /(47r) , replaced by 



(0) 



P a b(/x/,Mr) = a s (fi r )P^ b 



(0) i„ !*f 



I'r 



+«s 3 (/V 2 ) P, 



(2) 
ab 



WoP£ 



(25) 



ln^r 

Mr 



+/3o 2 P a ( b 0) ln 



2^f 



Mr 



where (3k represent the MS expansion coefficients 
of the /3-function of QCD |44I45I46147| . 

In Fig. the relative scale uncertainties of the 
average results is plotted, which is conventionally 
estimated by 



A/ = 



_ max [f(x, ~ min [f(x, /4)] 



2 | average [f(x,(4)} 



(26) 



where the scale $ varies p^. G [^^,4^^]. 

The top row of Fig. [S] shows the renormaliza- 
tion scale uncertainty of the NLO and NNLO pre- 
dictions for the scale derivative of q^ s , i = ±,v, 
as obtained from the quantity Aq^s defined in 
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NNLO : 
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5 10 4 


10 3 10 2 10 1 




X 



Figure 6. The relative scale uncertainty A/ (I26[) 
for all non-singlet and singlet cases. 



Eq. H'26|) ■ The bottom row displays the same, 
but for the singlet-quark distributions q s and the 
gluon distribution g. For all quark distribution, 
these uncertainty estimates amount to 2% or less 
at x > 10 -2 (for the gluon distribution 1% at 
x > 3 ■ 10~ 4 ), an improvement by more than a 
factor of three with respect to the corresponding 
NLO results. 

In general, for x > 10~ 3 the perturbative 
expansion for the scale derivatives / = 
d]nf(x ) ^if)/dln^f, f = q l ns , q s , g appears to 
be very well convergent and suggests a residual 
higher-order uncertainty of about 1% or less at 
a s < 0.2. Consequently the perturbative evolu- 
tion can be safely extended to considerably larger 
values of a s , hence lower scales, in this range of x. 
Larger corrections have to be expected at small x, 
but the results of the small- :r resummation alone 
will not help here. Further progress at small x 
would require at least a four-loop generalization 
of the fixed- N calculations |9I1UI11| and of the 



x-space approximations |43| linking them to the 
small- a; limits. In addition, one should also keep 
in mind that at fourth order the new colour struc- 
ture d abc d a bc/n c also will appear in singlet split- 
ting functions. 

4. CONCLUSION 

We have calculated the complete third-order 
contributions to the splitting functions governing 
the evolution of unpolarizcd parton distribution 
in perturbative QCD. 

The calculation is performed in Mcllin-iV space 
and follows the previous fixed-integer N compu- 
tations |9I1(JI11| . The extension to the complete 
analytical TV-dependence is the crucial new fea- 
ture. It required the set-up of an elaborate reduc- 
tion scheme for the corresponding loop integrals, 
an improved understanding of the mathematics 
of harmonic sums, difference equations and har- 
monic polylogarithms, and finally the implemen- 
tation of corresponding tools, together with other 
new features [Hj, in the symbolic manipulation 
program FORM 

Furthermore, by keeping terms of order e° in 
dimensional regularization throughout the calcu- 
lation, we have also obtained the third-order co- 
efficient functions for the structure functions F2 
and Fl in electromagnetic and for F3 in charged- 
current deep-inelastic scattering j^HJ. Addition- 
ally, the present method can be used to general- 
ize our fixed- N three-loop calculation of the pho- 
ton structure 0^1 to all N and it should also be 
possible to obtain the polarized NNLO splitting 
functions in this manner. 

The results for the three-loop splitting func- 
tions have been presented in both Mellin-A^ and 
Bjorken-a; space in Refs. |1I2I3| and agree with all 
partial results available in the literature. We have 
investigated the numerical impact of the three- 
loop (NNLO) contributions on the evolution of 
the parton distributions. The perturbative ex- 
pansion appears to be very well convergent ex- 
cept for very small x and shows good stability 
under variation of the scales. Thus, with the re- 
sults presented, the precision of the perturbative 
predictions for parton evolution has been greatly 
improvement. 
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